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Abstract. Two different computational approaches to the magnetorotational instability (MRI) are 
pursued: the shearing box approach which is suited for local simulations and the embedding box 
approach whereby a Taylor Couette flow is embedded in a box so that numerical problems with 
the coordinate singularity are avoided. New shearing box simulations are presented and differences 
between regular and hyperviscosity are discussed. Preliminary simulations of spherical nonlinear 
Taylor Couette flow in an embedding box are presented and the effects of an axial field on the 
background flow are studied. 



INTRODUCTION 

It is now generally accepted that the magnetorotational instability (MRI) is the main 
agent driving turbulence in accretion discs |1]. This instability provides a good expla- 
nation for the turbulence in accretion discs which might otherwise be hydrodynami- 
cally stable 10] • The historical developments of the understanding of this instability have 
appropriately been reviewed elsewhere in this book. Here we focus mainly on results 
within the shearing sheet approximation, which is the relevant extension of a periodic 
geometry to a shearing environment. Like with periodic "boundary" conditions, there is 
actually no boundary, and every point in the domain is equivalent to any other point. This 
avoids boundary layers which is extremely useful for simulations, but it is also astro- 
physically more relevant because we do not want explicit boundaries in a local domain 
that is supposed to represent a subvolume within the global disc. In the following we 
present a discussion of the nonaxisymmetric MRI in the shearing sheet approximation 
using the Rayleigh quotient to define in a convenient way an "instantaneous" growth rate 
of otherwise only transient growth. In this paper we also calculate the nonlinear evolu- 
tion of the MRI, reviewing both old simulations using subgrid scale modeling and new 
direct simulations where the ordinary viscosity and diffusion operators are used with 
uniform coefficients. Finally we turn to the issue of spherical Couette flow and present 
preliminary simulations in order to address recent experimental results suggestive of 
MRI in liquid sodium. 



AXISYMMETRIC VS NONAXISYMMETRIC MRI 



In order to obtain the dispersion relation for the MRI with a vertical magnetic field 
it suffices to consider the linearized pressureless one-dimensional momentum equation 
and the linearized induction equation for the two components in the rotational {x,y) 
plane, i.e. 

Ux - 2D.Uy = vaK, (1) 
Uy + {2- q)Q.Ux = VAb'y, (2) 

bx = VA^x, (3) 

by + qQ.b'x = VAby, (4) 

where dots and primes denote derivatives with respect to t and z, respectively, u and b are 
linearized velocity and magnetic field, and q denotes the steepness of the rotation law 
with Q. ~ R where <5r = 3/2 for purely keplerian discs. In our local cartesian model, y 
is the streamwise direction and x is the cross-stream direction. 

One could have retained the component of the momentum equation where the 
pressure gradient would enter, together with the continuity equation, but these two 
equations decouple from those considered already. The resulting new modes are the 
fast magnetosonic waves, which are of no particular interest in the present context. 

Assuming that the solution is proportional to e^'^^^^^', we obtain the dispersion relation 
jmi] in the form 

(0^ - 20)2 ^vlk^ + (2 - q)^^] + v\k^ {vlk^ - 2qQ}) = 0. (5) 

This is a bi-quadratic equation with altogether 4 solutions, corresponding to two different 
branches where ft) can have either sign on each of them. The two branches are 

(ol = vlk^ + {2-q)a'^± n^^vlk^ + (2 - q)^a'^. (6) 

The upper branch corresponds to Alfven waves and the lower branch corresponds to 
slow magnetosonic waves. The fast magnetosonic waves have been eliminated by using 
the pressureless momentum equation. In the range < v\k^ < 3 the slow magnetosonic 
waves can become unstable, i.e. co^ < 0, corresponding to an exponentially growing 
solution with growth rate Imft)_. The maximum growth rate is reached when v\k^ = 
and the corresponding value of ft)^ is then Imft)_ = ^Q.; see Fig.[ll 
In the nonaxi symmetric case the situation is not quite as simple, because now d /dy ^ 
0, and this leads to the occurrence of a new term with a non-constant coefficient. This 
comes from the U ■ V advection terms in all four equations, because U = Sxy. There are 
two different ways of dealing with this problem. One possibility is to abandon the as- 
sumption of harmonic solutions in the .x; direction, i.e. one assumes u = u(x)e^^y^'^^^~^^^^' . 
This approach has been taken by [4], for example. Another possibility is to assume 
shearing-periodic boundary conditions. The way then to get rid of the x dependence in 
the Sxd I dy term is by assuming k^ to depend on time, because then the d /dt term pulls 
down a factor k^x that can be arranged such as to cancel the Sxky term that results from 
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FIGURE 1. Solutions of the dispersion relation for accretion disc parameters (q = 3/2, thick solid 
lines), compared with the dispersion relation for a rigidly rotating disc (q = 0, dotted lines) and an almost 
nonrotating disc (q = and Q. = O.Olilo, dashed lines). Here, Q.q = GM/R? is the keplerian value at radius 
R. The upper set of curves denotes Alfven waves and the lower set slow magnetosonic waves. Note that 
in the completely nonrotating case (£2 = 0) the Alfven and slow magnetosonic waves are degenerate. 



the Sxd/dy term. (Here, as before, the dots denote a derivative with respect to t.) This 
requires k^x + Sxky = 0, i.e. 

k^{t)=k^Q-Stky, (7) 
where k^ is some initial values of k^. Using the ansatz 

q{x,y,z.,t) =AR&q{t) cx^[ik^{t)x + \kyy + \kzz\ , (8) 

where the hats denote the shearing sheet expansion, A is an amplitude factor, and 

^ Ay ^ ^ rrf 

q = [Ux, Uy, Uz, bx, by, bz, A)^ is the state vector with A = Cslnp,b = B/ y/jAQpo, the partial 
differential equation 

i^ + iSx^=Lq, (9) 

with non-constant coefficients, turns into an ordinary differential equation with constant 
coefficients 

i^=Lq. (10) 
at 

Without going into the details of the derivation we simply state the governing matrix, L, 
for the case with uniform density po, a uniform magnetic field in the y direction, yBo, 



and Alfven speed va = Bq/ ^/jMiPo and an isothermal equation of state, 
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(11) 



where we have used the abbreviations 2£L^ = 2£L-\- S, kf" = fc,VA, k'^ = kiC^. In the 
stratified case with a uniform background field, the system of governing equations has 
variable coefficients. Therefore we only consider here the case of constant density. Note 
that L is hermitian if 5 = 0. This property ensures that all eigenvalue are stable in that 
case. 

The set of equations (flOb were already investigated by using numerical integra- 
tion. They looked at the evolution q and found a transient behavior that depends on 
initial conditions. |Q, El studied the nonaxisymmetric stability in the presence of den- 
sity stratification which also gives rise to the Parker instability. A convenient method to 
calculate the growth rates of the MRI is in terms of the Rayleigh quotient (see [|2l]) 
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where {a\b) = ^a*bi defines a scalar product. To ensure solenoidality of the magnetic 
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field, we calculate bx for the initial perturbation from [6] 

bx{0)=~{kyby+k,b,)/kM. 



(13) 



The results for the keplerian case, S/Q. = —3/2, are shown in Figs|2land|3l In agreement 
with earlier work ll^ llflll . the maximum growth rate agrees with the Oort A- value lITlll . 
which is —5/2, or for keplerian rotation; see also Fig.[T] 

Thus, we see that in the nonaxisymmetric case the MRI shows great similarity with 
the axisymmetric counterpart. However, in the present case of the shearing sheet approx- 
imation it is strictly speaking only a transient. This becomes more obvious when looking 
at the evolution of see Fig. H] where we see an increase over about three orders 

of magnitude. In this figure we also show the resulting evolution of the root-mean- square 
magnetic field from a direct simulation of the shearing sheet equations, where we also 
adopt an isothermal gas with constant sound speed. The continuity equation is written 
in terms of the logarithmic density, A = Inp, 



(14) 
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FIGURE 2. Evolution of the imaginary part of co (a) and the norm of q (b) for vj^ky/£l = 1. Note that 
transient ampHfication is only possible during the time interval when \vj^k^^{t)\/Q.lQ. 
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FIGURE 3. Growth rates as a function of Alfven speed, in units of where k^^ denotes ky in 

the nonaxisymmetric case and k^ in the axisymmetric one. Symbols mark the results obtained from the 
Rayleigh quotient method for the nonaxisymmetric instability with ky = [0.8, 1.2]. For vj^ky/Q. > 1.7 our 
technique fails to yield reliable values of (O and the noisy oscillations seen in Fig.|2lbecome dominant. 
The solid lines indicates the result for the axisymmetric MRI (Eq.|6j. 
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and the induction equation is solved in terms of the magnetic vector potential A, where 
fi = V X A, and 

^A 

— = MxB + |^2oA^f+7]V2A, (15) 

r\ being the magnetic diffusivity and being the background rotation at the reference 
radius (the derivation of the shear term in this form is given in II2I1 I The momentum 
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FIGURE 4. Transient amplification of the magnetic field by the nonaxisymmetric magnetorotational 
instability. The solid line shows the result from the Rayleigh quotient method while the broken lines give 
the result from direct three-dimensional simulations with zero viscosity and zero resistivity. The square 
root of the Rayleigh quotient has been scaled by a factor Cg = 3.5 x 10^^ to make it overlap with Bnns 
curve. A resolution of only 32^ meshpoints is completely insufficient to resolve even the beginning of 
the instability. At least 256^ meshpoints are required to resolve the maximum. After fif > 17 even the 
simulation with 256^ meshpoints becomes under-resolved. 



equation is solved in the form 

_ = -M . Vm - c2 VA + — - + Fvisc + /, (16) 

where jQJt = d /dt -\-uf^ ^ and uf^ = —{3/2)Q.qx is the velocity in the y-direction 
due to the shear flow. Furthermore 7 = V x B//io is the current density, B is the magnetic 
field, fiQ is the vacuum permeability and F^isc 

is the viscous force. 

The initial condition for the 3-dimensional direct simulation is obtained by evolving 
linearized shearing sheet equations for ky = \ and = 10 to the point where fcx(^o) = ~5. 
[The size of the domain is (2;r)^.] For definitiveness, we reproduce here the numerical 
values in equation (lU): 

/-0.3108-0.0366i\ /+0.0683 - 0.5813i\ 

u= -0.4610-0.05421 , b=[ -0.0393 + 0.33401 , (17) 
V -0.0883 - 0.01041/ V +0.0347 - 0.29501/ 

and the logarithmic density is given by A = 0.042 — 0.36471. The amplitude is chosen to 
be A = 10^"^. Figure |5] shows images of on the periphery of the simulation domain. 




FIGURE 5. Images of the vertical component of the magnetic field, B^, for different values of t, where 
the field at f = corresponds to that given by equation ( I17> . 

The simulations have been carried out using the PENCIL CODE^ which is a high-order 
finite-difference code (sixth order in space and third order in time) for solving the 
compressible hydromagnetic equations. 

The way how this transient amplification can lead to sustained growth is through 
mode coupling, which is not considered in the present analysis. Relevant mode couplings 
could come about either through nonuniformities in the cross-stream or x direction and 
boundary conditions, or through nonlinearities. Li the shearing sheet approximation only 
the latter seems a viable possibility, and this is probably the mechanism through which 
the early shearing sheet simulations (e.g. Il4lll2|]) produced sustained turbulence. 



LOCAL DISC SIMULATIONS AND DYNAMOS 

The magneto-rotational instability is believed to be of great importance in connection 
with accretion discs. Here, gas spirals gradually into the center. This is possible mainly 
because of magnetic stresses, in particular those resulting from the small scale fields, 
brb^. They tap potential energy which gets converted into kinetic and magnetic energies, 
and eventually into heat which gets radiated away. This resistively produced radiation 
can be so big that it can explain the extreme radiation from quasars that are a hundred 
times more luminous than ordinary galaxies (e.g. lUsIl '). 

The standard theory of disc accretion is quite straightforward provided the discs can 
be treated as geometrically thin (e.g. lH^ The details of the magnetic stress are 
assumed not to matter. Based on dimensional reasoning the sum of the Reynolds and 
Maxwell stresses should scale like 

[u;iT0-b^/{^opo)] ^ ass^^'^H, (18) 

where ass is the magnetic contribution to the dimensionless Shakura-Sunyaev parameter 
its'], E = Jq pdz = 2poH, Q. = ^/gm]W is the keplerian velocity at radius R, and H is 
the scale height. In practice the Maxwell stress is always larger than the Reynolds stress. 
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but both tend to add to the stress, so there is no cancelation. Much of the work on the 
MRI in discs has focussed on determining the time averaged value of ass- In Fig.l^we 
show the evolution of ass for Run A of [ 19]. 

Here, much of the time variability comes from the fact that this simulation develops a 
large scale field, B, that is oscillatory. The dependence ass(B) can roughly be described 
via 

ass(B)^ag) + a^f)BV5^, (19) 

where ^ 0.002, a^g^ ^ 0.06, and Bq = /io(p)(Cs). In the original simulations of 
|l2] this time dependence of the large scale field was associated with an aO. dynamo. 
This large scale field also showed migration away from the equatorial plane. Various 
approaches have subsequently confirmed the initial result that a is negative in the 
northern hemisphere. This is quite unusual and is probably related to the swirl that 
comes from the combined action of magnetic buoyancy and the Coriolis force acting 
on the field-aligned flows driven by the B ■ VB tension force lE(jl . 

As expected, the dynamo al^ha requires the presence of both rotation and stratifica- 
tion. The simulations of iEH |22] without imposed field did show dynamo action, but 
they assumed periodicity in the vertical direction and there was no net stratification and 
hence no aO. dynamo-type behavior. Subsequent simulations of lEsIl had strong strat- 
ification and open boundary conditions. The vertical extent of the box was also much 
larger than in the previous simulation (|z| < 5H). Nevertheless, no large scale field was 
observed. At the moment we can only speculate about the possible origin of this dis- 
crepancy. One possibility is the fact that they used a numerical resistivity that allowed 
very little magnetic helicity to be dissipated. Subsequent work in the context of helically 
forced turbulence showed that a large scale dynamo effect involving the so-called a^ 
mechanism (for homogeneous helical turbulence, but no shear) can only saturate on a 
resistive time scale [24]. Whether or not this is indeed the reason for the absence of a 
large scale field in the simulations of lEsIl remains open. 



HIGH RESOLUTION DIRECT SIMULATIONS 



Several groups have shown independently that the MRI leads to sustained fully three- 
dimensional turbulence in the presence of an imposed (vertical or toroidal) magnetic 
field lfl3L Similar simulations have also shown that there is the possibility that 
the MRI coupled to the dynamo instability can lead to a doubly-positive feedback 
whereby the dynamo produces the magnetic field for the MRI, which in turn produces 
the turbulence for the dynamo EEUHQ. 

In all the simulation results published so far, a nonuniform effective viscosity or 
some other equivalent numerical procedure has been adopted in order to maximize the 
Reynolds number in regions where the flow is quiescent and to reduce it as much as 
necessary in regions where the flow shows strong spatial variations. These statements 
also apply to the effective magnetic diffusivity. Although these procedures are known 
to provide reasonably safe approximations to many types of turbulent flows [26], there 
are examples in the context of helical turbulence where departures from the ordinary 
viscosity and magnetic diffusion operators can lead to major differences compared with 
the correct solution [27] . 

In accretion discs the main source of heating is viscous and Joule heating. In the 
absence of cooling, this would lead to secular heating of the model. In order to avoid 
this, it is customary to add a volume cooling lfl9l] . In the present model we adopt instead 
an isothermal gas with constant sound speed. We thus solve equations (fT4ll - ([T6b using 
as initial conditions a simple magnetic field configuration given by 

= AQCOs{kix) cos{kiy) cos{kiz) : (20) 

where Aq = 0.2 has been chosen for all our runs. The initial field strength is chosen 
rather high so that the wavelength of the most unstable mode, va/^, is large compared 
with the mesh spacing, 5x. Alternatively, especially for high resolution runs, we use a 
snapshot from a lower resolution run and remesh it using interpolation. 

Unless noted otherwise, we present the results by measuring time in units of 
length in units of k^^, and density in units of the initial density, po. The orbital period 
is Trot = 2Tt/Q., which is sometimes also used when the time is given as t/T^ox- We 
have carried out simulations at three different resolutions with three different viscosities, 
keeping the magnetic Prandtl number Pr^ = v/r] always equal to unity; see Table[l] 

In the run with 128^ meshpoints we have been able to run for more than three hundred 
orbital times; see Fig.|7]for a plot showing the evolution of kinetic and magnetic energies. 
After the first 20 orbits the magnetic energy (per unit volume) decays rapidly and comes 
then to a halt at around 0.03 (the value quoted in Table[T]). 

Both kinetic and magnetic energies vary approximately periodically in intensity, but 
then at ? ~ lOOTrot the kinetic energy of the turbulence drops drastically by two orders of 
magnitude. The magnetic field also drops at first, but starts then to grow approximately 
quadratically in time, corresponding to a linear increase in the rms field strength. This 
increase is readily explained by the presence of a residual cross-stream rms magnetic 
field. Since in the simulation the magnetic vector potential is periodic (or shearing- 
periodic in the x direction), there is no net flux, so the most slowly decaying mode 
has the wavenumber k/ki = 1. This mode decays at a rate r^kf ^ 0.0015. This is 



TABLE 1. Mean magnetic and kinetic energies (per unit volume and in units of Eq = po^l^/k^) for 
runs at different resolution and different viscosities (for unit magnetic Prandtl number, v = Tj in units of 
n/fcj). The nondimensional stress is here given as ass = [{uxUy) — (bxby) /{nopo)]kl/Q?. For runs with 
regular viscosity, the magnetic Reynolds number is defined as = Mrms/(j?^i)- The duration of the run is 
indicated to give some idea about the statistical significance of the data. In all runs the normalized sound 
speed is Cski/Q. = 5. 
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FIGURE 7. Evolution of kinetic and magnetic energies in a run with 128 meshpoints. V = Tj = 
1.5 X lQ-^Q./kj, where ki = I = In/L and L = 2;r is the side length of the simulation box with uniform 
aspect ratio. 



indeed consistent with the data (see Fig. [51). In Fig.[5]we also demonstrate that the time 
integrated cross-stream rms magnetic field provides a rough estimate of the resulting 
toroidal magnetic field by winding up the poloidal field. Eventually, the field strength 
reaches an amplitude that is large enough for the MRI to set in (t ~ 200rrot). This leads 
to renewed turbulence for some 50 orbits, but then the magnetic and kinetic energies 
have dropped so much that the MRI becomes unable to sustain the turbulence (see Fig. [7] 
at t ^ 200rrot). 

To make contact with earlier work we now use hyperviscosity, i.e. the differential 
operator is replaced by V^. In numerical turbulence, hyperviscosity has previously 
been found to give rise to an artificially strong bottleneck effect, i.e. a shallower spectrum 
than k^^/^ near the dissipative subrange. However, more recent work has shown that 
the width of the bottleneck is always around one order of magnitude both for direct 
and hyperviscous simulations, and that the actual inertial range is not affected by the 




FIGURE 8. Comparison of (B?) 1/2 with {Bfj^l"^ (left panel) and / IQ.{bI) ^l^-At (shown as dash-dotted 
line) with {ByYl^ (right panel). 128^ meshpoints. V = 77 = 1.5 x lQ'^Q./k\. 




FIGURE 9. Evolution of kinetic and magnetic energies in runs using hyperviscosity and hyperresistiv- 
ity with 32^ meshpoints, V3 = 173 = 3.5 x \Q-^Q.Ik\ (left), and 16^ meshpoints, V3 = 173 = 7.5 x 10"^il//fcf 
(right). 



presence of hyperviscosity f2^. 

Using this as general reassurance of the usefulness of turbulence simulations with 
hyperviscosity, we now proceed using this technique in the present case of shearing 
sheet accretion disc turbulence. In the following, however, we shall only consider a 
rather small number of meshpoints, in which case there is no inertial range anyway. 
Therefore, such hyperviscous simulations cannot be considered as an approximation 
to the full equations, but they should really only be regarded as an illustrative model 
resembling features of a high resolution model. At a resolution of just 32^ we have in 
this way been able to run for 300 hundred orbits and found a steady level of turbulence 
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FIGURE 10. Kinetic and magnetic energy spectra for the run with 512 meshpoints, v = Tj = 2.5 x 
IQ-'^Q./k] (left) and 256^ meshpoints, v = Tj = 5 x \Q-'^Q./k\ (right). 

(left hand panel of Fig.|51). The absence of dramatic variations of the turbulence intensity 
suggests that such variations were an artifact of still too low resolution in the direct 
simulations presented above. This is supported by another lower resolution run (only 
16^ meshpoints) with hyperviscosity and hyperdiffusivity shown in the right hand panel 
of Fig. 121 We see that the behavior is similar to what we see in Fig.|7J suggesting that the 
reason for the disappearance of the MRI is really just that the Reynolds number dropped 
below a certain critical value. 

In the rest of this section we consider the properties of the energy spectra for the 
direct simulations with high resolution of 256-^ and 512-^ meshpoints. It turns out that 
the magnetic energy spectrum increases with wavenumber like k^/^ . The kinetic energy 
scales approximately like k^^^^ . For neither of the two spectra do we have a reasonable 
explanation. Obviously, the larger the resolution, the harder it becomes to run for a long 
time. Our run with 256^ meshpoints has run for nearly 120 orbits, but this may not be 
enough to exclude conclusively a subsequent decay. 

The simulations presented above may help to determine a critical value of the mag- 
netic Reynolds number, i?m"'\ above which dynamo-generated turbulence becomes pos- 
sible. If the turbulence in the simulation with 256-^ meshpoints remains indeed sustained, 

(crit) 

we expect R)^ to be between 100 and 200. However, the values of ocss iiot con- 
verged in any of the simulations in that its value increases with resolution, and it is also 
much larger in the simulations with hyperviscosity than with ordinary viscosity. 

SPHERICAL COUETTE FLOW (PRELIMINARY) 

In order to better understand the experimental verification of the MRI in spherical 
Taylor-Couette flow by the Maryland group i28ll . it is desirable to perform simulations 
that are related to the setup used in the laboratory. We know already that cylindrical 



FIGURE 11. Toroidal magnetic field component displayed on the periphery of the computational 
domain (color coded). 512-' meshpoints. v = Tj = 2.5 x IQ^'^Q./k^. 



Taylor-Couette flow can weU be modeled by embedding the flow in a cartesian box 
1122] • We also know that simulations of the geodynamo (convection-driven dynamos in a 
spherical shell) can successfully be modeled by embedding the spherical shell in a box 
|30]. It is therefore natural to attempt modeling of the Maryland MRI experiment using 
the same setup. In Fig. [121 we present some preliminary results of spherical Couette flow 
embedded in a cartesian box with an imposed axial magnetic field of given strength. 
We work in SI units and chose parameters that are closely related to the experimental 
setup. A sphere of radius 7? = 0.15 m is embedded in a box of size (2 X 0.18m)3. The 
sphere spins with an angular frequency of Q. = 300 s^^ Outside the sphere the velocity 
is damped to zero at rate T^|^p = 200 s^^ The density is initially uniform and equal 

to p = 930kg m^-^. For reasons of computational simplicity our simulation is weakly 
compressible, and we use an isothermal equation of state (ratio of specific heats is 7= 1) 
with sound speed Cs = 50ms^. For the present simulations no attempt is made to use 
realistic values of the kinematic viscosity and the magnetic diffusivity; both are chosen 
equally big with v = 7] = 2 x 10^-^m^s^^ The field strength is varied between 500 G 
and 2kG. 

It turns out that for weak and strong fields the magnetic field is completely axisym- 
metric. For an intermediate field strength of 5o = 1 kO both the flow and the magnetic 
field become nonaxisymmetric. 

As anticipated by R. HoUerbach (private communication), the main reason this sim- 
ulation produces nonaxisymmetric flows is the modification of the background flow in 
such a way that it becomes Rayleigh unstable, i.e. q = — dlnQ/dlnr > 2. On the other 
hand, this does not seem to be the case in the Maryland experiment where the flow 
remains close to keplerian, i.e. q ^ 1.5. 
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FIGURE 12. Images of the vertical component of the magnetic field, Bv, for three different values of 
the imposed field, Bo- The images shown on the two side walls of the box go actually through the axis 
(x — y — 0), the top plane goes through z ~ 0.09m and the bottom plane goes through z = 0. 



The main outcome of the magnetorotational instability is by now well known and its 
importance for astrophysical discs is obvious. Many details, however, remain far from 
clear due to the fact that many of the numerical shearing sheet simulations invoke arti- 
ficial viscosity which is hard to quantify. Thus, we have no clear idea about the critical 
value of the magnetic Reynolds number required for the combined magnetorotational 
and dynamo instabilities. The present work suggests that the number is around 200, and 
certainly above 100, but obviously a more quantitative study would be desirable. Fur- 
thermore, details of the nonlinear mode interactions are not well understood. This would 
be very useful for obtaining a clearer picture of how in the shearing sheet approximation 
the nonaxisymmetric MRI leads to sustained growth. The Pencil Code is well suited 
for such studies, and it is publicly available, so it may be hoped that the remaining gaps 
in our understanding will soon be filled by new members of the scientific community. 

Regarding spherical Couette flow, the present approach using embedded box simula- 
tions may not be optimal, and dedicated codes for spherical geometry may be advan- 
tageous. However, any independent verification using different methods is still often 
useful. Returning to the astrophysical context, it is important to assess the relative im- 
portance of dynamo-generated fields relative to those generated on a more global scale in 
the rest of the disc. This question can only be appropriately in a fully global simulation. 



CONCLUSIONS 



Indeed, significant progress has recently been made in this direction 
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